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The main goal of numerical relativity is the long time simulation of highly nonlinear spacetimes that cannot be 
treated by perturbation theory. There are three elements to achieving this: 

(1) Analytic issues, such as well-posedness, constraints, boundary conditions, linear stability, gauge conditions and 
singularity avoidance. 

(2) Computational issues, such as evolution and boundary algorithms, numerical stability, consistency, spacetime 
I discretization and numerical dissipation. 

(3) Physical issues, such as simulation of the desired global spacetime, extraction of the radiation from an isolated 
O ' system, the proper choice of initial data, long timescale evolutions tracking many orbits of an inspiraling binary. 

The correct treatment of the physical issues introduces severe global problems. For instance, long term simulations 
are needed to flush out the spurious gravitational radiation contained in the initial data for the gravitational field 
of a binary black hole so that a physically relevant waveform can be extracted. Furthermore, extraction of the 
waveform requires a compactified grid extending to null infinity, or some alternative approximation based upon an 
I I outer boundary in the radiation zone. 
. At present, the major impasses to achieving such global simulations are of an analytical/computational nature. 
£N| ' We present here some examples of how analytic insight can lend useful guidance for the improvement of numerical 
ON . approaches. 
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II. WAVES 

The prime physical objective is to compute the gravitational waves emanating from a compact source. We begin 



by introducing some underlying mathematical and computational problems in terms of two examples of scalar wave 
propagation. Both of these examples are chosen because they have a direct analogue in general relativity and illustrate 
computational problems that arise because of exponentially growing modes in an analytic problem which is well-posed. 



A. Unbounded exponential growth 

First, consider a nonlinear wave propagating in Minkowski space according to 

r) af3 d a d $ - ij7 Q/3 (9 Q $)(^$) = = $J7 Q/3 9 Q ^log$. (2.1) 

Although this nonlinear equation arises from a linear wave equation for log<I>, it is a remarkably accurate model for 
some of the problems that occur in numerical relativity. In order to simplify the problem, we first impose periodic 
boundary conditions so that the evolution takes place in a 3-torus T 3 , i.e. on a boundary free manifold. For $ > 
the Cauchy problem for this system is well-posed. Furthermore, the linear superposition log $1 + log $2 of solutions 
to the linear wave equation correspond to the solution $i < i>2 of the nonlinear problem. 
A nonsingular solution of this system is the wave 

® = l + F{t-z), (2.2) 

where F > —1. Suppose we try to simulate this solution numerically. If numerical error excites an exponentially 
growing mode of this system then noise in this mode will eventually dominate the wave being simulated. For the 
linear wave equation there are no such exponential modes. But this nonlinear system admits the solutions 

$ A = e At (l + F(t - z)), (2.3) 
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for arbitrary A. 

Thus, although we have a well-posed initial value problem whose principle part is the Minkowski wave operator, the 
simulation of a simple traveling wave is complicated by the existence of neighboring solutions which grow exponentially 
in time. Numerical error will excite these exponentially growing modes and eventually dominate the traveling wave 
we are attempting to simulate. You might ask: Why not choose log $ as the evolution variable? This would clearly 
solve all numerical problems. As we have already said, we have chosen this example because it arises in numerical 
relativity where there is no analogous way to take the logarithm of the metric. But there are indirect ways to model 
the derivative of a logarithm, analogous to grouping terms according to and rewriting the nonlinear wave 

equation (2.1) as 



This indeed works for the scalar field, as illustrated in Fig. 2. (See the discussion of finite difference methods in 
Sec. III.) We will come back to the gravitational version of this problem but first we give another example which 
illustrates a similar complication with the initial-boundary value problem. 

We base this example on the initial-boundary value problem (IB VP) for this nonlinear scalar field in the region 
—L < z < L obtained by opening up the 3-torus to T 2 x R. We consider the simulation of a traveling wave wave 
packet $ = 1 + F(t — z) > with the Neumann boundary condition 8 z $\ z= ±l = d z F\ z= ±[ J . The wave packet gets 
the correct Neumann boundary data for it to enter the boundary at z = —L, propagate across the grid and exit the 
boundary at z = +L. 

In the process, numerical noise will be generated. There are solutions of the system of the form 



for arbitrary e > 0. Normally, if a scalar field admitted such solutions we could infer that the corresponding Cauchy 
problem was ill-posed by arguing (following Hadamard) that the solution $ = has vanishing Cauchy data and that 
the neighboring solutions <& e have unbounded size for any t > 0. However, the above Cauchy problem is well-posed 
because the initial data must satisfy <f> > 0. 

After the wave packet has crossed the grid, the remnant numerical noise gets homogeneous Neumann data 
d z &\ z= ±L = 0. Thus it is reflected off the boundaries and trapped in the simulation domain where it can grow 
exponentially. The noise is generated while the wave packet is traveling across the grid. The short wavelength modes 
can be controlled by introducing numerical dissipation. But the long wavelength modes cannot be damped without 
the danger of damping the signal. Just as in the case of periodic boundary conditions, numerical error can excite 
exponential modes that destroy the accuracy of a simulation in the case of Neumann boundary conditions. 

You might ask: Why use Neumann boundary conditions? The Sommerfeld boundary condition (<9 t ±d z )$\ z= ±L = 
does not admit such modes and moreover it propagates numerical noise off the grid. The answer to that question has 
to wait until we have discussed the constraint equations of general relativity. 

The problem with using Neumann boundary conditions is not of analytic origin. The problem is of a numerical 
nature. Whereas the signal gets the correct inhomogeneous boundary data to propagate it off the evolution domain, 
the noise gets the left-over homogeneous data and gets trapped in exponentially growing modes. 

The lesson here is that it is preferable to use Sommerfeld type boundary conditions, not for physical or mathematical 
advantage but for numerical advantage. A Sommerfeld boundary condition doesn't solve all the problems. Even 
though a homogeneous Sommerfeld condition carries energy out of the evolution domain, it does not in general 
give the physically correct outer boundary condition for an isolated system. For a nonlinear system such as general 
relativity, one would need an inhomogeneous Sommerfeld condition whose boundary data could only be determined 
by matching to an exterior solution. But numerically a Sommerfeld condition has the great advantage of allowing the 
noise to escape through the boundary. Unfortunately, in present numerical relativity codes, Sommerfeld boundary 
conditions are inconsistent with the constraints, which we will discuss later. 



Another mechanism by which a reflecting boundary condition can introduce exponential modes is the repetitive blue 
shifting off moving boundaries. This can even happen for a linear wave propagating between two plane boundaries 
in Minkowski space. The boundaries can effectively play ping-pong with a wave packet by arranging the boundary 
motion to be always toward the packet during reflection. 



V ^d a (-d f3 ^>) = 0. 



(2.4) 




(2.5) 



B. Moving boundaries 
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Let x a = (i,x,y,z) be inertia coordinates, with the reflecting boundaries in the (x,y) plane. Under reflection, 
functional dependence of a wave packet traveling in the positive z-direction changes according to 

$(£-z)^$(e 2a (t + z)), (2.6) 

where the speed of the reflecting wall is 

w^tanha. (2.7) 

After many reflections the energy in the wave grows by a factor e 4aT , where T is measured in units of the crossing-time 
between reflections. 

It is instructive to reinterpret this experiment from a numerical relativity viewpoint where the spatial coordinates 
of the boundaries have fixed grid values. For that purpose, we consider the well-posedness of the IBVP for the linear 
wave equation 

g^dadp* = 0. (2.8) 

in a general background spacetime with non-constant metric g a 0, Again let the evolution domain be the region 
-L < z < L. 

Most of the mathematical literature on well-posedness of the IBVP is based upon symmetric hyperbolic systems 
in first derivative form [1-4]. We achieve this for the wave equation by introducing auxiliary variables u = ($, d a $>). 
Standard results then imply a well posed IBVP for a homogeneous boundary condition of the matrix form Mu = 
provided that 

• the resulting energy flux normal to the boundary has the dissipative property 

^"(u) > (2.9) 



• M has maximal rank consistent with this dissipative property and 

• M is independent of u. 

In the present case, the energy flux is determined by the energy momentum tensor for the scalar field. The flux 
normal to the boundary at z = +L is 

T n = -n a {d t <$>)d a <5> (2.10) 

where n a — g ZOL /\Jg zz is the unit normal to the boundary. 

The dissipative condition can be satisfied in a variety of ways. The choice 

<9 t $ = (2.11) 

leads to a homogeneous Dirichlet boundary condition; and the choice 

n a d a <5> = 0. (2.12) 

leads to a homogeneous Neumann boundary condition. Homogeneous Dirichlet and Neumann boundary conditions 
are limiting cases for which T n = 0, i.e. there is no energy flux across the boundary and signals are reflected. Between 
these limiting cases, there is a range of homogeneous boundary conditions with T n > of the form n a d a <& + Pdt& = 0, 
where P > 0. Of particular interest is the Sommerfeld-like case where the derivative lies in an outgoing characteristic 
direction. 

The IBVP for the scalar wave equation is well-posed for any of these maximally dissipative boundary conditions. 
Furthermore, by consideration of the symmetric hyperbolic equation satisfied by u — q(x a ), where q has explicitly 
prescribed space-time dependence, the well-posedness of the IBVP with the homogeneous boundary condition Mu = 
can be extended to the inhomogeneous form M(u — q(t, x, y)) = 0, with freely assigned boundary data q. By using 
such boundary data, a Neumann or Dirichlet boundary condition can be used to model a wave which is completely 
transmitted across the boundary with no reflected component, at least at the analytic level. 

Note the important feature that the free boundary data for the scalar field consist of one function of three variables 
in contrast to two functions for free Cauchy data. As we shall see, this is the major complication in formulating a 
constraint preserving boundary condition for a well-posed IBVP in general relativity. 
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Now consider the simulation of a linear scalar wave in the flat background metric which results from the transfor- 
mation t = t, x = x, y = y, z = z + A sin from inertial coordinates x a . In these x a coordinates, the boundaries at 
z = ±L are oscillating relative to the inertial frame, as indicated by the "shift" g zt — — Aucosojt. For our simulation, 
we prescribe initial data $o = d t $o = and either the appropriate Neumann or Dirichlet boundary data q~L{t) and 
q+h{t) for a wave packet which enters the boundary at z = —L, travels across the domain and leaves the boundary 
at z — +L. A second order accurate code would simulate this signal with 0(A 2 ) error in the grid displacement A. 
Thus <& = 0(A 2 ) after the packet has traversed the domain. However, this remnant error gets homogeneous boundary 
data. Although, as just discussed, the normal energy flux associated with homogeneous Dirichlet or Neumann data 
vanishes in the rest frame of the boundary, in the x a inertial frame the boundary is moving and the noise can be 
repeatedly blue shifted, resulting in an exponential increase of energy. 

One way to eliminate this problem would be to deal with coordinate systems in which the shift vanishes, at least at 
the boundary. However, this would rule out many promising strategies for dealing with binary black holes, e.g. the 
use of co-rotating coordinates or of generalized Kerr-Schild coordinates. But especially in a nonlinear problem such 
as general relativity, the excitation of exponential modes can rapidly destroy code performance. 

III. GENERAL RELATIVITY: HARMONIC EVOLUTION 

The previous examples of scalar waves show that even if the underlying analytic problem is well-posed and even 
if the numerical simulation converges to the analytic solution, the existence of exponentially growing modes in the 
analytic system can effectively invalidate long term code performance. In general relativity, coordinate freedom is 
a further complication which can introduce rapidly growing modes that are an artifact of gauge pathologies. In 
order to illustrate computational problems that are not a trivial consequence of gauge, we consider the harmonic 
formulation of Einstein's equations. Although no coordinates can guarantee complete avoidance of gauge problems, 
harmonic coordinates have several advantages for investigating the interface between numerical and analytic problems 
in general relativity: 

• Small number of variables 

• Small number of constraints (4 harmonic conditions) 

• Einstein's equations reduce to quasilinear wave equations 

• Well-posed Cauchy problem [5] 

• Symmetric hyperbolic formulation [6] 

• Global asymptotically flat solutions for weak Cauchy data [7] 

• Well-posed homogeneous IBVP [8] 

A numerical code for evolving Einstein's equations, the Abigel code [8], has been based upon a generalized version 
of harmonic coordinates satisfying the curved space wave equation 

H a := ^—gU g x a = d^^—gg^d,x a ) = H a (x^g pry ), (3.1) 

where H a are harmonic source terms. These harmonic source terms do not affect any of the analytic results regarding 
well-posedness but, in principle, they allow any spacetime to be simulated in some generalized harmonic coordinate 
system. The harmonic reduced evolution equations are written in terms of the metric density -f^ u — ^/^gg^" whose 
ten components obey quasilinear wave equations 

l^dadp^ = S" v , (3.2) 

where S M " contains nonlinear first-derivative terms that do not enter the principle part. The harmonic conditions 
C M := — = are the constraints on this evolution system which are sufficient to ensure that Einstein's 
equations are satisfied. Except where noted, we set = to simplify the discussion but all results generalize to 
include nonvanishing gauge source terms. For details concerning the formulation and its implementation see Ref. [8]. 
By virtue of the evolution equations, the harmonic constraints satisfy a homogeneous wave equation of the form 

7 Q/3 dadpC + A^ a d a C p + B^C P = 0. (3.3) 
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Thus, in the domain of dependence of the Cauchy problem, the solution C M = is implied by standard uniqueness 
theorems provided the system is initialized correctly. 

A well-posed evolution system is a necessary but not sufficient ingredient for building a reliable evolution code. 
Code performance can be best tested by simulating an exact solution and measuring an error norm. The error should 
converge to zero in the continuum limit as the grid spacing A shrinks to zero. In testing evolution codes it is desirable 
to first eliminate effects of boundary conditions by imposing periodicity in space, which is equivalent to carrying out 
the simulation on a 3-torus without boundary. A suite of toroidal testbeds for numerical relativity has been developed 
as part of the Apples with Apples project [9,10]. The convergence and stability of several codes [9,11,17], including 
the Abigel code [9], has been demonstrated using this test suite. 

One testbed is the Apples with Apples periodic gauge wave with metric 

ds 2 = <I>(-dt 2 +dz 2 ) + dx 2 + dy 2 , (3.4) 

where 

* = 1 + Asin(?^). (3.5) 
It is obtained from the Minkowski metric ds 2 = —di 2 + dx 2 + dy 2 + dz 2 by the harmonic coordinate transformation 




V = V- 




FIG. 1. Snap shots of a gauge wave simulation carried out with an early version of the Abigel code. The code is stable and 
the error converges to zero at second order in grid spacing A but after a few crossing times the error is too large to make the 
results useful. 
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Figure 1 shows some snap shots of a gauge wave simulation carried out with an early version of the Abigel code. 
The time dependence of the error shows exponential growth of the form 



£ ~ A 2 iP(t,z)e M . (3.7) 

Here ip is a well-behaved function which is almost identical in shape to the signal As a result, the error cannot 
be dissipated by standard techniques for dealing with short wavelength noise. The exponential growth originates in 
exact analogy with our example for the nonlinear scalar wave equation (2.1). In fact, the metric (3.4) is a flat solution 
of the Einstein equations in harmonic coordinates for any z) which satisfies the nonlinear wave equation (2.1). 
The general solution is <f> = e /(*- 2 0+s , (*+ z ). j n particular, there are exponentially growing harmonic gauge waves 

ds 2 = <5>x(-dt 2 +dz 2 ) + dx 2 + dy 2 , (3.8) 

where 

^ = e-(l + Asin(^ r ^)), (3.9) 

which lie arbitrarily close to the gauge wave being simulated. Thus numerical error inevitably excites exponential 
modes which eventually dominate the simulation of the gauge wave. The practicality of code performance depends 
on the timescale of this exponential growth. 

Although the Abigel code is stable, convergent and based upon a well-posed initial value problem, like any other code 
it is subject to the excitation of exponential modes in the underlying analytic system. Certain numerical techniques 
greatly improve its accuracy in simulating the gauge wave: 

• Expressing the equations into flux conservative form, an idea from computational fluid dynamics which was 
introduced into general relativity by the Palma group [16]. 

• Summation by parts, introduced into general relativity by the LSU group [17], which at the level of linearized 
equations leads to energy estimates for the semi-discrete system of ODE's in time which arise from spatial 
discretization. 

• Nonlinear multi-pole conservation, which suppresses the excitation of long wavelength exponential modes by 
grouping the troublesome nonlinear terms in a way that enforces global semi-discrete conservation laws (or 
approximate conservation laws). 

The semi-discrete multipole technique (being introduced here) provides an excellent example of how analytic insight 
into the source of a numerical problem can be used to design a remedy. As we will show, various combinations of 
these three techniques lead to dramatic improvement in gauge wave simulations. Other numerical methods based 
upon enforcing or damping the constraints are not crucial for the gauge wave problem but can be important for 
simulations of curved spacetimes. 

An essential ingredient in any code is the method used to approximate derivatives. The Abigel code treats the 
quasilinear wave equations (3.2) as first differential order in time and second order in space. This allows use of explicit 
finite difference methods to deal with the mixed space-time derivatives introduced by the "shift" term in the wave 
operator while avoiding the artificial constraints that would be introduced by full reduction to a first order system. 
On a grid with spacing A, the natural finite difference representation for the the first and second spatial derivatives 
are the centered approximations 

W )^)= ^ + A) 2 / ( *- A) (3.10) 

and 

31F {Z ) - DVFiz) = F( - + A) - 2 ^ ) + F( - A) . (3.11) 

These formulae were used in the simulation labeled TIGHT in Fig. 3. Although the code was tested to be stable 
and convergent with second order accuracy, the excitation of the exponentially growing mode of the analytic problem 
limits accurate simulations to about 10 crossing times on a reasonably sized grid. 

In order to exaggerate nonlinear effects, the simulations shown in Fig. 3 were carried out for a highly nonlinear 
gauge wave with amplitude A = .5, on a scale where the metric is singular for A = 1. (The standard Apples with 
Apples tests specify amplitudes of A = .01 and A = .1.) Problems with exponential modes do not appear for small 



G 



amplitudes simulations in the linear domain. One contributing factor to the exponential growth is that the tight 3- 
point stencil (3.11) for the second derivative does not lead to an exact finite difference representation of the integration 
by parts rule necessary to establish energy conservation, which is the main idea behind the summation by parts (SBP) 
method. But this is only part of the story since standard SBP techniques only apply to linear systems. 

It is instructive to examine how these ideas extend to the second derivative form of the nonlinear wave equation 
(2.1) which underlies the gauge wave problem. This will illustrate in a simple way how flux conservative equations, 
SBP and multipole conservation can combine to suppress excitation modes in the analytic problem. The model 
scalar problem is effective in isolating the difficulties underlying a full general relativistic code, in addition to allowing 
efficient computational experimentation. 

We carry out the analysis for waves traveling with periodic boundary conditions in one spatial dimension. The 
extension to three dimensions is straightforward but notationally more complicated. The theory regarding well- 
posedness of hyperbolic systems is based upon the principle part of the equations. For that reason, we first consider 
the linear wave equation 

<9 Q <9 Q $ = -<9 2 $ + <9 2 $ = 0. (3.12) 
The energy associated with this wave can be related to the conserved integral 

T = $ 2 ] = ^ (*i0"$2 - $2^*i)dV^ (3.13) 
by choosing <&i = <& and $2 = <9t$. For the case of periodic boundary conditions on the interval < z < L, 

I = J o ^(<9 t $) 2 -$<9 2 $^(fe. (3.14) 

The integration by parts by parts rule 

,-L / 

d z {<Pd z <S>) + (d z $) 2 + <Pd 2 z $)dz = 0, ( 3 - 15 ) 



applied to a periodic interval, then supplies the key step in using the wave equation to relate I to the positive definite 
energy 

1 = £ = £ (W) 2 + {d z ^dz. (3.16) 

In order to obtain a discrete version of the integration by parts identity (3.15), we introduce a uniform grid Zi, 
< i < N , with spacing A and represent 



/' 

Jo 



L N 

Fdz^ Aj2fi+i/2 (3.17) 



where 



_ F( Zi ) + F(z i+1 ) 
Ji+1/2 2 ' *■ ' 



In addition we represent derivatives at the midpoints by the centered approximation 



d z F{ Zl + A/2) - f i+1/2 = F{Zl+1 \ F{Zl) (3.19) 



so that periodic boundary conditions imply 



la o 
This ensures the semi-discrete monopole conservation law 



/ d z Fdz^Aj2.n +1/2 =F\^=0. (3.20) 

JO n 



dl <b $dz -> 0, (3.21) 
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which results from the flux conservative form of Eq. (3.12). Equation (3.21) controls growth of the spatial average of 
$ but not of its non-constant spatial Fourier components which measure its gradient. 

Energy estimates control the growth of the gradient of $. With the above definitions, it is straightforward to check 
that 

d z (FG) — Fd z G — Gd z F — > 0. (3.22) 

As a result, the semi-discrete version of the integral identity (3.15) is satisfied if the second derivative term is rep- 
resented as a product of first derivatives. For the linear wave equation this results in the semi-discrete conservation 
law 

d t £ -> 0. (3.23) 

In order to implement SBP in a code such as the Abigel code, which is second differential order in space with the 
fields represented by their values on grid points, and not on mid-points, the above results can be applied by treating 
the mid-points for even numbered grid points as the odd- numbered grid points, and vice versa. This results in the 
widened finite difference representation for the second derivative 

*«.) - D 2 F{z) = DDF(z) = ^ + 2A)-2FM + F(.-^ &M) 

as opposed to the tight stencil (3.11). Figure 2 shows the remarkable improvement in long term accuracy obtained 
in the simulation of a non-linear wave satisfying Eq. (2.1). (Numerical dissipation has been used to damp short 
wavelength instabilities triggered by the loose coupling between even and odd grid points.) The curves labeled 
TIGHT are obtained using the standard stencil (3.11). They show exponential growth on a scale of ss 10 grid crossing 
times. The curves labeled SBP are obtained using the stencil (3.24) consistent with SBP. This change of stencil 
suppresses growth of long wavelength exponential modes so that accurate simulations of rs 1000 crossing times are 
possible. 



TIGHT 

- SBP 

- SBP-MON 

- TIGHT-MON 




time (crossing times) 

FIG. 2. A comparison of the various evolution algorithms used to evolve the nonlinear wave equation (2.1). The tests are 
based on the sine wave solution (3.5), with amplitude A = .5, simulated on a grid of N = 100 points with a time-step of 
At — Az/4. The graph shows the l x norm of the error. 



Summation by parts only has general applicability to linear equations, although the technique extends in an ap- 
proximate sense to the nonlinear domain. Other approaches can also be successful for nonlinear problems, especially 
if the troublesome nonlinear terms can be identified. In the case of the nonlinear equation (2.1), these terms can be 
incorporated in the principle part by reformulating the equations in the flux conservative form 

d a (v a ^d A=O, (3.25) 
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with the subsequent reduction reduction to the first order in time system 



and 



d t <$> = $Q (3.26) 
d t Q = d z (^d z $). (3.27) 

Many choices of spatial discretization of this flux conservative system lead to an exact semi-discrete version of the 
monopole conservation law 

d t [ Q = 0. (3.28) 
Jo 

We consider the two choices 

<m>i.-KH, +1/2 -HH- 1/2 - (3 - 30) 

As a result of either of these discretizations, the initial data determine the conserved value of the monopole moment 
J Q L Qdz and the excitation of the exponential modulation (2.3) is thereby frozen out of the numerical evolution. In 
this way a tight 3-point stencil (3.30) can be used, as opposed to the wide 5-point stencil (3.29) (and the concomitant 
numerical dissipation) required by SBP for the second order system. The curve labeled TIGHT-MON in Fig. 2 
shows how long term accuracy is dramatically enhanced by this technique, without use of numerical dissipation. The 
curve labeled SBP-MON shows that, in this case, no additional improvement is gained when monopole conservation 
is combined with SBP. 

These numerical techniques introduced for the model scalar problem were formulated in a way that could be 
readily taken over to the gravitational case. Although the Einstein equations can neither be linearized by taking the 
"logarithm of the metric" nor written in a completely flux conservative form analogous to (3.25), there are various 
ways to group derivatives which decouple the Jacobean transformation that generates the exponential mode in the 
gauge wave metric (3.8). One example is the grouping 

g^d"g^ = (6?5l + 5°^) J-^$ A) (3.31) 

for which expression of the principle part of the Einstein equation in the form d p {g a ^d p g fl ij) leads to the semi-discrete 
conservation laws 

dt I g^&gtfdz -> (3.32) 
Jo 

for the gauge wave. The conserved quantities are comprised of multipoles of monopole (the spatial trace) and 
quadrupole (the trace-free part) type. 

The advantage of enforcing these conservation laws is exhibited in Fig. 3. Comparison of Figs. 2 and 3 shows 
that SBP and multipole conservation lead to almost identically beneficial results in simulating the gauge wave as 
in simulating the nonlinear scalar wave. The standard 3-point stencil (TIGHT) again excites exponentially growing 
error on the order of 10 crossing times but accurate evolutions for over 1000 crossing times are attained either with 
SPB or with a 3-point stencil embodying multipole conservation (TIGHT-MULT). 



9 



1.5 



TIGHT 

- SBP 
SBP-MULT 

- TIGHT-MULT 




200 400 600 800 1000 

time (crossing times) 

FIG. 3. A comparison of the various evolution algorithms used to evolve the harmonic Einstein equations. In these tests the 
code evolved flat spacetime in the gauge defined by Eq. (3.6), with amplitude A — .5. The size of the grid was N = 100, with 
a time-step of At = Az/4. The graph shows the £oo error norm of the g zz metric component. 




time (crossing times) 

FIG. 4. A comparison of the various evolution algorithms used to evolve the harmonic Einstein equations. In these tests the 
code evolved the gauge wave metric with shift defined in Eq. (3.33), with amplitude A — .5. The size of the grid was N = 100, 
with a time-step of At = Az/4. 



Whether objectionable gauge modes can be decoupled so effectively in a more general problem is an interesting 
question. However, all the above numerical techniques, which lead to excellent code performance for the gauge wave, 
are of a universal nature that can be adopted for the simulation of a general spacetime by a general code. Since most 
simulations contain a weak field region, such as the far field region outside a black hole, these techniques might in 
fact be necessary in order to avoid excitation of local versions of exponential Minkowski gauge modes. Figure 4 shows 



10 



how these methods extend to the challenging simulation of the gauge wave with shift 



ds 2 = -(1 - Asina)dt 2 + 2Asinadtdz + (1 + A sin a) dz 2 + dx 2 + dy 2 , 



(3.33) 



with a — w(t — z)/L. The simulation was carried out with periodic boundary conditions and amplitude A = .5, so 
that the grid has an effective velocity of half the speed of light. Again there are exponentially growing gauge waves, 



(for arbitrary A), which satisfy these boundary conditions. These will trigger a numerical instability unless their 
excitation is controlled by a conservation law on the semi-discrete system. Remarkable improvement in long term 
performance is achieved by implementing either SBP or the multipolc algorithm. 

These examples show what must be done, beyond having a stable, convergent code, in order to achieve accurate long 
term simulations. Exponential modes undoubtedly arise in a wide variety of systems with the examples presented 
here just the tip of the iceberg. Short wavelength modes arising from discretization error can be suppressed by 
numerical dissipation. The long wavelength modes exist in the analytic problem. This raises some key questions: Arc 
there geometric clues to identify the origin of such long wavelength exponential modes? What numerical or analytic 
techniques can be used to suppress them? 



Given an evolution code on the 3-torus which is based upon a well-posed Cauchy problem for Einstein's equations 
and which is free of all numerical problems, several things can go wrong in extending the evolution to include a 
boundary. On the analytic side, the imposition of the boundary condition can be ill-posed or it can lead to violation 
of the constraints or it can introduce exponentially growing modes. On the numerical side, the finite difference 
implementation of the boundary condition can be unstable or inaccurate. On the physical side, the correct boundary 
data representing radiation (or the absence of radiation) entering the system might not be known or it might not be 
possible to extract the waveform of the outgoing radiation. 

Here we examine the analytical and numerical issues for the harmonic IBVP. The reduced evolution system consists 
of the quasilinear wave equations (3.2). Our discussion for nonlinear scalar waves show that the IBVP for this system 
is well-posed for any maximally dissipative boundary conditions, e.g. Dirichlct, Sommcrfeld or Neumann. 

Next consider the harmonic constraints C^. They satisfy the homogeneous wave equation (3.3). Thus we can 
formulate a well-posed IBVP for the propagation of the constraints by imposing a maximally dissipative boundary 
condition. Then, given that the constraints and their time derivative are satisfied by the initial data and that the 
constraints have homogeneous boundary data, the uniqueness of the solution to the constraint propagation equations 
would imply that the constraints be satisfied in the domain of dependence of the IBVP. However, consistency between 
the boundary conditions for the evolution variables and the homogeneous boundary conditions for the constraints is 
not straightforward to arrange. 

For example, consider evolution in the domain z < with boundary at z = 0. In the tangential-normal 3+1 
decomposition x^ = (x a , z) intrinsic to the boundary, a homogeneous Dirichlet condition on the constraints takes the 
explicit form 



A naive attempt to satisfy these conditions by boundary data on the evolution variables would involve assigning 
both Dirichlet (tangential) and Neumann (normal) conditions to -f az , which would be an inconsistent boundary value 
problem. 

One way to impose consistent constraint preserving boundary conditions is based upon the well-posedness of the 
Cauchy problem. Consider smooth Cauchy data which is locally reflection symmetric with respect to the boundary 
at z = 0. Then in some neighborhood— L < z < L of the hypersurface z = the Cauchy problem is well-posed. On 
z = 0, the local reflection symmetry implies that the evolution equations satisfy 



ds\ = -(e At - Asina)dt 2 + 2,4 sin adtdz + (e xt + Asina)dz 2 + dx 2 + dy 2 



(3.34) 



IV. THE HARMONIC IBVP 



C z = d al za + d zl zz = 
C a = d bl ab + d zl az = 0. 



(4-1) 
(4.2) 



7 za = 
d zl zz = 
d zl ab = 



(4.3) 



and that the constraints satisfy 
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C z = 
d z C a = 0. 



(4.4) 



It is straightforward (although algebraically complicated) to show for the harmonic IBVP that the combination 
of Dirichlct and Neumann boundary conditions (4.3) implies that the constraints satisfy the homogeneous boundary 
conditions (4.4). Thus (4.3) provide homogeneous constraint preserving boundary conditions for a well-posed harmonic 
IBVP. 

Well-posedness of the IBVP extends to the case of "small" boundary data, of the form M(u — q(x a )) = discussed 
in Sec. Ill, in the sense that the prescribed data q is linearized off a solution with homogeneous boundary data. 
However, the available mathematical theorems do not guarantee well-posedness for finite boundary data. We describe 
below the major issues regarding constraint preserving inhomogeneous boundary conditions for the harmonic IBVP. 
For further details, see Ref. [8]. 

Part of the inhomogeneous boundary data which generalize (4.3) are associated with the gauge freedom corre- 
sponding to a boundary version of the "shift". By a harmonic coordinate transformation it is alway possible to 
set 

Y- a = q a (x h )Y- Z (4.5) 

at the boundary, where q a is freely prescribed data. The unit normal to the boundary then defines the normal 
derivative 

&n ■■= j^"^ = d * + ( 4 - 6 ) 

entering the Neumann boundary data, q zz = d n -y zz and q ab — d n j ab , which complete the inhomogeneous version of 
(4.3). 

The boundary data q = (q a , q zz , q ab ) can be freely prescribed in a well-posed IBVP for the reduced evolution system 
but they must be restricted to satisfy (4.4) in order to ensure that the constraints are satisfied. The condition C z = 
requires 

q zz = ~d a q a j ZZ - (4.7) 

When the boundary shift — q a is nonzero, the second condition in (4.4) must be restated in the form d n C a = 0, 
because the derivative d z is no longer in the normal direction to the boundary. This condition is a restriction on the 
data q ab , which are closely related to the extrinsic curvature K ab of the boundary. It requires that 

(,22 

V^hD b (K b - 5 b a K) + VF^afcC" - y —C b d a q b = 0, (4.8) 

where h a b and D a are the metric and connection intrinsic to the boundary. This equation can be recast as a symmetric 
hyperbolic boundary system which determines the 6 pieces of Neumann data q ab in terms of 3 free functions, the free 
(gauge) data q a and the boundary values of r y zz , j ab and d z j za . Solutions of reduced equations with this boundary data 
necessarily satisfy the constraints. Unfortunately, the appearance of the quantities 7 ZZ , r y ab and d z j za complicates 
the well-posedness of the constrained IBVP since these quantities cannot be freely specified but must be determined 
in the course of the evolution. 

Formally, the constraint preserving boundary data have the functional dependence q = q(u, x a ), which involves 
evolution variables u whose boundary values cannot be freely prescribed. This complication has its geometric origins 
in the fact that the boundary data (gauge quantities and extrinsic curvature) do not include the intrinsic metric, as 
in the case of Cauchy data. Because of the dependence of the constraint preserving boundary data on u, the available 
theorems regarding well-posedness only apply to perturbations of homogeneous data, where the background values of 
u can be explicitly determined. 

These constraint preserving boundary conditions have been implemented in the Abigel code. Test simulations of 
the IBVP for the shifted gauge wave (3.33) were carried out by opening one face of the 3-torus to form a T 2 x [0, 1] 
manifold with boundary. Figure 5 shows the results reported for an early version of the code [8] . The graphs indicate 
stability and and convergence but there is also a growing error which eventually leads to a nonlinear instability. One 
underlying cause of this error growth is the continuous blue shifting off the moving boundaries, as discussed in Sec. II. 
However, these tests were carried out before semi-discrete conservation laws were incorporated into the evolution 
algorithm so that a better understanding of the error must await future test runs. 
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FIG. 5. The norm of the finite-difference error 7" = jana ~ Inum, rescaled by a factor of 1/A 2 , for a gauge-wave. The 
tests were carried out with an early version of the Abigel code before semi-discrete conservation laws were incorporated. The 
upper two (mostly overlapping) curves demonstrate convergence to the analytic solution for a wave with amplitude A = 10 _1 
gridsizes 80 3 and 120 3 . We also plot |i?|oa, the £00 norm of {H 1 ) 2 + S^H^-Hi , to demonstrate that convergence of the 
harmonic constraints is enforced by the boundary conditions. The lower curve represents evolution of the same gauge-wave 
with A = 10 -3 for 300 crossing times, with gridsize 80 3 . 



V. SOMMERFELD ALTERNATIVES 



The examples presented here indicate a computational advantage in formulating boundary conditions in a manner 
such that numerical noise can propagate off the grid for the case of homogeneous boundary data. To date, there exists 
only one well-posed formulation of the IBVP for general relativity that allows this type of generalized Sommerfeld 
boundary condition. This is the Friedrich-Nagy formulation [13] based upon a formulation of Einstein's equations in 
which an orthonormal tetrad, the connection and the Weyl curvature are treated as evolution variables. The gauge 
freedom in the theory is adapted in a special way to the boundary so that boundary conditions need only be imposed 
on the curvature variables. The critical feature of the formalism is that the constraints propagate tangential to the 
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boundary. This allows the well-posedncss of the IBVP for the reduced evolution system to be extended to the fully 
constrained system. Unfortunately, this formulation has not yet been implemented as a numerical code, partially 
because of its analytic complexity and partially because it would require some infrastructure beyond that existing in 
most present codes. 

An important issue is whether this success of the Friedrich-Nagy system in handling a Sommerfeld boundary 
condition is limited to formulations that include the tetrad or the curvature among the basic evolution variables. In 
linearized gravitational theory, there is a simple variant of the harmonic formulation that has a well-posed IBVP, 
admits a Sommerfeld boundary condition and has been successfully implemented computationally [15]. The nonlinear 
counterpart consists of the evolution system 

-f f, d a df,i i *=S i 1, (5.1) 
H a := d tl ta + drf" = H a (x,j), (5.2) 

comprised of the wave equations (3.2) for the six spatial components jij and propagation equations for the time 
components -f ta . Alternatively, the propagation equations could be reformulated as 

d t H a = d t H a (5.3) 

in order to make the evolution system uniformly second differential order. Wcll-posedness of the nonlinear Cauchy 
problem does not follow in any direct way from standard theorems. An analysis of the principle part shows that 
this naive harmonic system is only weakly hyperbolic, which opens the door for lower derivative terms to produce 
instabilities [14]. 

It is instructive to investigate the performance of a code based upon this weakly hyperbolic harmonic system by 
using the Apples with Apples testbed. Figure 6 shows the results of the robust stability test, where a simulation in 
the linear regime is carried out with random (constraint violating) initial data. The results show an exponential rise 
in the violation of the Hamiltonian constraint at a rate that increases with grid resolution, which eventually leads 
to a code crash. This behavior is symptomatic of weakly hyperbolic systems and presages possible problems in the 
nonlinear domain. The simulation of a nonlinear gauge wave with shift, shown in Fig. 7, verifies such problems. 
These problems do not appear for the nonlinear gauge wave without shift, as the results shown in Fig. 8 indicate 
convergence. Also, as illustrated in Fig. 9, with the addition of numerical dissipation, the Hamiltonian constraint 
no longer grows exponentially in the robust stability test, although the constraint violation still increases with grid 
resolution, indicating failure of the test. Similar conclusions follow from the Apples with Apples Gowdy wave tests. 

These results show that a full battery of tests are necessary in order to establish reliable code performance. Other- 
wise, misleading information about code performance can result. As history has shown in the case of ADM evolution 
codes, weakly hyperbolic systems system cannot be expected to give reliable long term performance in the presence 
of strong fields, which makes them unsuitable for black hole simulations. 
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FIG. 6. The robust stability test for the weakly hyperbolic harmonic system. The £ x norm of the Hamiltonian constraint is 
plotted on a linear-logarithmic scale. All specifications are in accord with the Apples With Apples test. 
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FIG. 7. The nonlinear gauge wave with shift test for the weakly hyperbolic harmonic system. The code crashes in less than 
a crossing time. 
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FIG. 8. The nonlinear gauge wave without shift test for the weakly hyperbolic harmonic system, run in accord with the 
Apples With Apples specifications. The convergence of the error is deceptive of code reliability. 
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FIG. 9. The robust stability test for the weakly hyperbolic harmonic system with dissipation. As in Fig. 6, the Hamiltonian 
constraint is plotted on a linear-logarithmic scale. The dissipation now kills the exponential growth but the growth of constraint 
violation with resolution indicates that the code fails the test. 



The Friedrich-Nagy system and the weakly hyperbolic harmonic system represent two extremes of a dilemma facing 
code development in numerical relativity. On one hand, the Friedrich-Nagy system has all the desired analytic features 
but its complexity poses a barrier to code development. On the other hand, the weakly hyperbolic harmonic system 
is simple and easily implemented as an efficient code, but well-posedness is questionable. Should you try to fix these 
simple systems or should you bite the bullet and develop codes based upon formulations where a well-posed nonlinear 
IBVP has been fully established? To date, the effort in numerical relativity has been weighted heavily toward the 
simpler formulations. It is timely that some attention be given to investigating whether the Friedrich-Nagy system 
can be converted into a workable code. A useful starting point would be a linearized version evolving on T 2 x R, 
where the complications of the equations and the boundary gauge would greatly simplify and would perhaps lead to 
a better understanding of the essential elements of the approach. Most of the effort in the field can be expected to 
remain a compromise between these extremes, e.g. the strongly hyperbolic harmonic system for which a Sommerfeld 
boundary condition is not constraint preserving. In all such endeavors, a close working combination of analytic and 
numerical insight can offer valuable guidance. 
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